      SUBROUTINE REDUCE (N,MA,A,MB,B,LOW,IGH,CSCALE,WK)
CSE
C
CPS       PURPOSE:
C
C             REDUCES, IF POSSIBLE, THE ORDER OF THE GENERALIZED
C             EIGENVALUE PROBLEM A*X = (LAMBDA)*B*X BY PERMUTING
C             THE ROWS AND COLUMNS OF A AND B SO THAT THEY EACH HAVE
C             THE FORM
C                        U  X  Y
C                        0  C  Z
C                        0  0  R
C
C             WHERE U AND R ARE UPPER TRIANGULAR AND C,X,Y,AND Z ARE
C             ARBITRARY. THUS, THE ISOLATED EIGENVALUES CORRESPONDING
C             TO THE TRIANGULAR MATRICES ARE OBTAINED BY A DIVISION,
C             LEAVING ONLY EIGENVALUES CORRESPONDING TO THE CENTER
C             MATRICES TO BE COMPUTED.
C
CPE
CAS          ***  ARGUMENT LIST  ***
C
C  INPUT PARAMETERS:
C
C  N - THE ORDER OF THE MATRICES A AND B;
C
C  MA - THE ROW (FIRST) DIMENSION OF THE ARRAY A AS SPECIFIED
C       IN THE CALLING PROGRAM;
C
C  A - A REAL TWO-DIMENSIONAL ARRAY WITH ROW DIMENSION MA AND COLUMN
C      DIMENSION AT LEAST N CONTAINING THE MATRIX A AS THE PROBLEM
C      IS DEFINED IN THE FUNCTION STATEMENT ABOVE;
C
C  MB - THE ROW (FIRST) DIMENSION OF THE ARRAY B AS SPECIFIED IN
C       THE CALLING PROGRAM;
C
C  B - A REAL TWO-DIMENSIONAL ARRAY WITH ROW DIMENSION MB AND COLUMN
C      DIMENSION AT LEAST N CONTAINING THE MATRIX B AS THE PROBLEM
C      IS DEFINED IN THE FUNCTION STATEMENT ABOVE.
C
C  OUTPUT PARAMETERS:
C
C  A - CONTAINS THE PERMUTED A MATRIX;
C
C  B - CONTAINS THE PERMUTED B MATRIX;
C
C  LOW - THE INTEGER SPECIFYING THE BEGINNING INDEX OF THE SUBMATRICES
C        OF A AND B CONTAINING THE NON-ISOLATED EIGENVALUES;
C
C  IGH - THE INTEGER SPECIFYING THE ENDING INDEX OF THE SUBMATRICES
C        OF A AND B CONTAINING THE NON-ISOLATED EIGENVALUES. IF
C        IGH=1 (THEN LOW MUST BE 1 ALSO), THE PERMUTED A AND B
C        MATRICES ARE UPPER TRIANGULAR;
C
C  CSCALE - THE REAL ARRAY OF DIMENSION AT LEAST N CONTAINING THE
C           REQUIRED COLUMN PERMUTATIONS IN ITS FIRST LOW-1 AND ITS
C           IGH+1 THROUGH N LOCATIONS;
C
C  WK - THE REAL ARRAY OF DIMENSION AT LEAST N CONTAINING THE
C       REQUIRED ROW PERMUTATIONS IN ITS FIRST LOW-1 AND ITS
C       IGH+1 THROUGH N LOCATIONS.
C
CAE
C  REQUIRED SUBPROGRAMS - NONE
C
C  REQUIRED FUNCTIONS - NONE
C
C  REFERENCE - R.C. WARD, BALANCING THE GENERALIZED EIGENVALUE
C              PROBLEM, TECHNICAL REPORT ORNL/CSD-47, OAK RIDGE,
C              TENNESSEE, 1979.
C
C              [SAME TITLE] SIAM J. SCI. STAT. COMPUT., VOL 2, 
C                           NO 2, JUNE 1981, PP 141-152.
C
C ****
C
      DIMENSION A(MA,N),B(MB,N),CSCALE(N),WK(N)
      DOUBLE PRECISION A,B,CSCALE,WK,F
      K = 1
      L = N
      GO TO 20
C ****
C     FIND ROW WITH ONE NONZERO IN COLUMNS 1 THROUGH L
C ****
   10 L = LM1
      IF (L .NE. 1) GO TO 20
      WK(1) = 1
      CSCALE(1) = 1
      GO TO 200
C
   20 LM1 = L-1
      DO 70 II=1,L
         I = L+1-II
         DO 30 J=1,LM1
            JP1 = J+1
            IF (A(I,J) .NE. 0. .OR. B(I,J) .NE. 0.) GO TO 40
   30    CONTINUE
         J = L
         GO TO 60
   40    DO 50 J=JP1,L
            IF (A(I,J) .NE. 0. .OR. B(I,J) .NE. 0.) GO TO 70
   50    CONTINUE
         J = JP1-1
   60    M = L
         IFLOW = 1
         GO TO 150
   70 CONTINUE
      GO TO 90
C ****
C     FIND COLUMN WITH ONE NONZERO IN ROWS K THROUGH N
C ****
   80 K = K+1
C
   90 DO 140 J=K,L
         DO 100 I=K,LM1
            IP1 = I+1
            IF (A(I,J) .NE. 0. .OR. B(I,J) .NE. 0.) GO TO 110
  100    CONTINUE
         I = L
         GO TO 130
  110    DO 120 I=IP1,L
            IF (A(I,J) .NE. 0. .OR. B(I,J) .NE. 0.) GO TO 140
  120    CONTINUE
         I = IP1-1
  130    M = K
         IFLOW = 2
         GO TO 150
  140 CONTINUE
      GO TO 200
C ****
C     PERMUTE ROWS M AND I
C ****
  150 WK(M) = I
      IF (I .EQ. M) GO TO 170
      DO 160 IS=K,N
         F = A(I,IS)
         A(I,IS) = A(M,IS)
         A(M,IS) = F
         F = B(I,IS)
         B(I,IS) = B(M,IS)
         B(M,IS) = F
  160 CONTINUE
C ****
C     PERMUTE COLUMNS M AND J
C ****
  170 CSCALE(M) = J
      IF (J .EQ. M) GO TO 190
      DO 180 IS=1,L
         F = A(IS,J)
         A(IS,J) = A(IS,M)
         A(IS,M) = F
         F = B(IS,J)
         B(IS,J) = B(IS,M)
         B(IS,M) = F
  180 CONTINUE
C
  190 GO TO (10,80), IFLOW
C ****
C     EXIT REDUCE
C ****
  200 LOW = K
      IGH = L
      RETURN
      END
